Forum for Electromagnetic Research Methods and Application Technologies (FERMAT) 


Convolution Quadrature and the Time Domain 
Integral Equations of Electromagnetics 


Jielin Li’, Peter Monk*, and Daniel S. Weile” 
‘Department of Electrical and Computer Engineering, University of Delaware 


"Department of Mathematical Sciences, University of Delaware 
(Email: weile@udel.edu) 


Abstract—The use of convolution quadrature (CQ) 
approaches for the discretization of time domain integral 
equations (TDIEs) is described in full. Using an 
operational calculus approach, CQ methods render the 
continuous TDIE convolution discrete through a 
mapping from the Laplace domain to the Z domain. This 
process simplifies the computation of the spatial 
integrations needed for the integral equation 
discretization, as the shadow region endemic to temporal 
Galerkin discretization is eschewed. The underlying 
frequency domain nature of CQ also eases its use for 
dispersive kernels. Numerical results will demonstrate 
the technique. 


Index Terms—Time domain integral equations (TDIE), 
Convolution Quadrature (CQ), marching on in time 
(MOT). 


I. INTRODUCTION 


Time domain integral equations (TDIEs) have 
long been the béte noir of computational 
electromagnetics research. Historically unstable, 
TDIEs have no commercial implementation to this 
day. Nonetheless, the past two decades have seen 
remarkable advances in TDIE implementation, 
application and stabilization, and today TDIEs are 
not uncommon subjects of research. 


The earliest methods for the discretization of 
time in TDIEs relied almost exclusively on 
temporal Galerkin (TG) methods, in which a single 
basis function is chosen and shifted repeatedly in 
time to represent the current on a given spatial 
basis function [1-4]. Unsurprisingly, this approach 
led to a concentration on the choice of basis 
function as a panacea for stability concerns, and a 
flurry of works suggesting candidate temporal 
basis functions on the basis of experimental 
evidence followed [5-8]. 


The folly of this quest was demonstrated through 
the work of Toufic Abboud and Jean-Claude 
Nédélec and other mathematicians, who solved the 
stability problem with careful spatial integration 
[9, 10]. When computing the electromagnetic 
interaction between spatial basis functions, this 
approach used a combination of analytic 
integration and a careful computation of the 
intersection between basis function support and 
illuminated regions. Previous methods for kernel 
element computation, regardless of the choice of 
basis function, used Gaussian (or other fixed point, 
fixed weight) integration rules defined over 
regions unconnected with the patch illumination. 
They thus yielded inaccurate results when sharp 
boundaries between light and shadow invalidated 
the assumptions underlying the derivation of the 
integration rule. 


The gradual comprehension of the source of the 
trouble lead to a burst of activity seeking to 
compute the elements correctly by smoothing the 
illuminated region/shadow’ region boundary. 
Thus, in the early part of the century, Weile et al. 
described the use of bandlimited interpolatory 
functions that made the shadow region boundary 
essentially continuous and differentiable to all 
orders [8]. Because these basis functions were 
necessarily noncausal, a complicated scheme for 
extrapolation and time marching were used to 
arrive at a stable scheme. 


Another way of avoiding the shadow region 
involves employing full-domain basis functions in 
time; the so-called marching on in degree (MOD) 
technique removes the shadows this way [11]. 
Finally, this type of smoothing approximation was 
brought to its logical conclusion by Pray et al., 
who expanded prosaic basis functions into series 
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of smooth functions that can be integrated easily, 
but still have causal, compact temporal support 
[12]. 

In this paper, we review an unrelated method of 
temporal discretization known as convolution 
quadrature (CQ) [13-16]. In the context of TDIE 
solution, CQ refers to the computation of the 
necessary convolution using a refined form of 
Heaviside’s operational calculus. In modern 
terms, the method consists in beginning with a 
Laplace domain equation, and discretizing it by 
replacing the Laplace domain variable s with a 
rational or matrix function of the discrete 
transform variable z to achieve a discretization in 
the frequency domain. This gives rise to 
expressions that have several computational 
advantages, and are provably stable in most 
instances of scientific import. In this paper, we 
describe the application of CQ to electromagnetics 
problems and give several examples. We begin 
with an overview of the method itself. 


II. CONVOLUTION QUADRATURE 


In this section, we discuss the convolution 
quadrature procedure itself. Because it involves 
mathematical manipulations that may not be 
familiar to practitioners of computational 
electromagnetics, we begin with a section of 
mathematical preliminaries. 


A. Mathematical preliminaries 


To begin to discuss the CQ approach to the 
discretization of integral equations, we review the 
Laplace and Z transforms, and a few fine points of 
matrix algebra. The Laplace transform of a 
function of time f(t) is a function 


f(s) = | Fedt, (1) 


where the integral converges. The lower limit of 
the integral is taken to be O so that any 
singularities located at the origin may be included. 
Moreover, if the function f(t) does not grow 
faster than exponentially, that is, if 


lim f(e-** = 0, (2) 


for some real a, then the integral converges for 
Re(s) < a. In particular, if the integral converges 
for all s with negative real parts (1.e., for all values 
in the left half plane), the function represented 1s 
bounded as t > œo [17]. 


For our purposes, the most important feature of 
this representation of a function is that quality 
most beloved of sophomore engineering students 
subjected to courses in differential equations: The 
Laplace transform converts differentiation to 
multiplication. Specifically, if g(t) = f'(t), then 


G(s) = sf (s) — f (07). (3) 
Indeed, this formula motivates most Laplace 
domain approaches, and even the transform 
nomenclature: by “transforming” our problem into 
a different “domain” we simplify theoretical 
manipulations by converting calculus into mere 
algebra. This observation becomes even more 
intuitive and powerful for physical applications, 
where it is assumed that f is continuous and 
strictly causal so that f (07) = 0. 


Less well known in_ the computational 
electromagnetics community is the Z transform, 
which functions as a sort of Laplace transform for 
discrete sequences [18]. Given a sequence 
fn, for n = 0,..., 0, its Z transform is given by 


f=) faz. (4) 


(Note that we assume all sequences are causal, 
that is, we define fa = 0 for all n < 0.) This form 
is immediately recognized as a Laurent series with 
coefficients taken from the series f, itself. 
Because of this, the coefficients can be recovered 
using Cauchy’s integral formula 


1 F =4 
f= zo ora O 
C 


where the contour C encircles the origin inside the 
region of convergence of the transform [18]. (if 
the sequence described is bounded for large values 
of n, this includes the entire z-plane outside the 
unit circle. The unit circle itself is included if 
fi =O(n%) for n>o and a >0.) This 
formula will give rise to very efficient methods for 
the computation of the CQ kernel. 
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Operational formulas like Eq. (3) also hold for 
the Z transform. In particular, it should be clear 
that defining the delayed sequence gn = fnim 
where m is a positive integer, results in the 
formula 


giz) =z "f(Z). (6) 
(Of course, this form of the formula requires the 
causality assumption mentioned above. A similar 
formula can also be written for sequences that are 
advanced, if care is taken to remove the initial 
values of the sequence that would otherwise move 
to negative indices.) 


Finally, we will have occasion to compute 
functions of (square) matrix argument [19]. This 
is not so daunting as it as it at first seems; 
consider, for instance, a monomial function 
f(A) = kA". This function is easily computed: 
after all, A can simply be multiplied by itself n 
times, and the result scaled by k. The method for 
computing an arbritrary function of A follows 
immediately: if the Taylor series for f(x) is 


CO 


fx) =) ax”, 0) 
n=0 
then 
f(A) = > anA”, (8) 
n=0 


and there is no ambiguity in this last formula 
(provided it converges). The required 
computation can be accelerated, however, if the 
eigendecomposition of A is computed. Suppose 
we write the matrix as 


A = XAX~! (9) 


where X is a vector of eigenvectors, and A is a 
diagonal matrix of eigenvalues. (We ignore for 
simplicity and computational realism the 
possibility that A has a geometrically repeated 
eigenvalue.) Substituting Eq. (9) into Eq. (8) 
shows that 


f(A) = Xf(A)X™. (10) 


Computing a function of a diagonal matrix is quite 
easy, since it merely involves computing the 
function of the diagonal elements as can be 
inferred from Eq. (8). 


With these preliminaries firmly in hand, we can 
now begin to describe the CQ method itself. 


B. Multistep Based CQ 


CQ discretization of integral equations begins 
with spatial discretization of the equation in the 
Laplace domain. Whether the spatial 
discretization is accomplished with the method of 
moments or the Nystrom method, and regardless 
of basis function selection, integration formula, or 
the equation being discretized, this process results 
in a matrix equation of the form [20-23] 


Z(s)I(s) = V(s). (11) 
Assuming N unknowns are used in the process, 
Z(s) is an NXN matrix representing the integral 
equation operator, V(s) is an N-vector related to 
the electromagnetic excitation, and I(s) is an N- 
vector of unknowns. 


Eq. (3) shows that multiplication of a function 
by s in the Laplace domain is tantamount to 
differentiating it in the time domain. Eq. (6) shows 
that multiplication by z~* in the Z domain is 
tantamount to a delay in the discrete time domain. 
Therefore, one way to effect a temporal 
discretization would be to substitute a finite 
difference approximation in z~* for s. 


A little reflection reveals that a more general 
approach is necessary; after all, the finite 
difference approximation to the derivative may 
also be assigned to different places in the interval. 
Given a time step At, and writing 


fn = f (nt) (12) 
we may decide, for instance, that the 
approximation 

n+1 ` fı n 
E A 13 
f N (13) 


approximates f'(nAt + At), f'(nAt + *At), or 
f'(nAt). This choice can be reflected in our 
discretization by taking a more general point of 
view that examines differential equations rather 
than merely derivatives. 

Thus, consider the differential equation 


dy 
= 14 
it = SY ( ) 
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Here “s” is a parameter, but this equation still 
relates the idea that multiplication by s is 
equivalent to differentiation [13]. (In particular, 
“s” is the ratio between the function’s derivative 
and its value.) One method of solving this 
equation numerically is called a linear multistep 
method, in which the unknown y is sampled at 
integer multiples of a time step At (i.e. y, = 
y(nAt)), and approximated iteratively using a 
formula of the form [24] 
k k 


2 AiYn+i-k = SAt » DY iG (15) 
i=0 


i=0 
While this formula may look unfamiliar in this 
form, it merely represents the most general way of 
solving Eq. (14) on a discretized grid of size At. 
For instance, if k = 1, a, = 1, a7 = —1, 641 = 0, 
and po = 1, the formula assigns a first-order finite 
difference approximation of the derivative to the 
rear-most point in the interval; that is, it uses Eq. 
(13) to represent the derivative at nAt, rather than 
any other point in the interval. In the numerical 
differential equations literature, this method is 
known as “forward Euler.” “Backward Euler,” or 
BE, obtains if k = 1, a, = 1, a7 = —1, 6 = 1, 
and fo = 0; this amounts to the same derivative 
approximation being assigned to (n + 1)At. The 
trapezoidal method approximates the derivative in 
the center of the interval, and so corresponds to 
k=1, a,=1, a =-1, Bp =fho= 5. For 
reasons that will become clear later, the method 
with the best combination of accuracy and stability 
properties for our purposes is the BDF2 (second 
order backward difference formula) method, with 


coefficients k =2, a,= Z 


2 
Bo = 1, and fp, = Bo = 0. 
Computing the Z transform of Eq. (15) in view 
of Eq. (6) suggests the substitution [13] 
SEn = 
At yar B,z'-* At 
This equation can be substituted into Eq. (11) and 


inverse Z transformed yielding the set of matrix- 
vector convolution equations 


= aL 
Ai = —2,Qo —_ >) 


(16) 


n-1 
Zoln = Vy — X A A (17) 
m=0 


This process can be illustrated analytically and 
in pictures. Consider, for simplicity, the 1x1 
Laplace domain matrix kernel 


Z(s) =e **, (18) 


where Tt > 0 is a positive constant. (This is, of 
course, merely the Laplace transform of the time 
domain kernel 6(t—T)). If this “matrix” is 
discretized with the backward Euler method, 
Q(z) = 1 — z™t, we obtain 


T T 
Z(z) = exp (- ac) exp (z) (19) 
By Taylor expansion, we see immediately 
lyr” 
— TANF ELSEN PENRE -n 
Bi a) a T Z (29) 
n=0 
so that 
1 tty” 
SEENI 21 
n is T eh 


These figures demonstrate that the outcome of CQ 
discretization 1s a discrete approximation of the 
delay corrupted by numerical dispersion. In 
integral equation implementations, this dispersion 
introduces error, but also makes the spatial 
integrals easier to compute. 


While the above description of the CQ method is 
straightforward, it leaves open the question of the 
stability of the method. Assuming that the original 
system Eq. (11) is stable, all singularities a of Z(s) 
will have the property that Re(a) <0, ie. all 
singularities of the system matrix will lie in the left 
half plane (LHP). (On the other hand, a pole in the 
right half plane would correspond to an 
exponential growth.) By the same token, the Z 
transform of a causal and stable system kernel has 
all its singularities inside the unit circle; that 1s, all 
singularities a of Z(z) have the property that 
la| < 1. Therefore, any mapping z = Q71+(sAt) 
that maps the LHP inside the unit circle will 
preserve the stability of Eq. (11) upon 
discretization [13]. 


Of course, the mapping z = Q~*(sAt) is merely 
the inverse of Eq. (16), derived from the original 
multistep method of Eq. (15). A multistep method 
with the property that its RHP is mapped into the 
unit circle is known as A-stable in the mathematics 
literature since it is stable independent of step size 
when used to solve differential equations [24]. 
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BE, BDF2, and the trapezoidal method are all A- 
stable, but the (forward) Euler method is not. 
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(b) 
Fig. 1. The response Zn of Eq. (21) with At = 1 for different T. 


While this observation implies that A-stability 1s 
necessary for an accurate and stable discretization 
of the integral equation, it is not quite sufficient. 
Some A-stable methods, when applied to stiff 
ordinary differential equation problems, have a 
distressing tendency to oscillate as time 
progresses. (In this context a stiff differential 
equation 1s one involving processes that evolve on 
very different time scales. Many numerical 
schemes encounter stability or efficiency 
difficulties solving stiff equations because the time 
steps employed must resolve fast processes while 
not unduly burdening the simulation of slow 
processes. ) 


For integral equation applications, this tendency 
appears as a kernel discretization that never 
decays. A multistep method is called L-stable if, 
when applied to Eq. (14) as Re(s) > —o, it yields 
a discretization that decays and vanishes at late 
time steps, rather than oscillating. BE and BDF2 
are L-stable, but the trapezoidal rule is not [24]. 


C. Kernel Computation 


Another important property of the kernel 
discretization method described in the previous 
subsection is that, if need be, it can be computed 
numerically with surprising alacrity and precision. 
Specifically, if the Laplace domain form of the 
kernel is known, the Z transform inversion 
described in Eq. (5) can be computed accurately 
and efficiently using the humble trapezoidal rule 
[13]. 


Suppose that an element of the Laplace domain 
Z matrix is called (s). Inserting Eq. (16) into Eq. 
(5), and choosing a circle of radius p < 1 for the 
contour C yields the inversion formula 


z 21T ja 
in = | 4 o ) eJ9"dQ, (22) 


0 


This computation will be most well conditioned if 
there are no singularities on the unit circle and the 
contour radius is chosen as unity. If this is not 
possible, the radius may be chosen as p = 1— ô 
for 6 « 1. In any case, because the integrand is 
periodic, this integral can be estimated with 
exponentially increasing accuracy as a function of 
M, the number of points used in its estimation. 

This yields the approximation 

z m 
p” © |e (pe?"™ ) jon 
in ap Dl ae 23) 
m=0 


This equation, in turn, is recognized immediately 
as the discrete Fourier transform of the sequence 


: m 
pr folo] og 
Bm = M’ At 
for m=0,..,M-—1 [18]. Therefore, this 
computation is very inexpensive—it can be 


Forum for Electromagnetic Research Methods and Application Technologies (FERMAT) 


performed with a fast Fourier transform (FFT), and 
converges exponentially once sampled adequately. 
Moreover, it is easily applied to very complicated 
Green’s functions that are not easily expressed in 
the time domain, such as are needed to compute 
interactions of dipoles through arbitrarily 
multilayered structures. Indeed, this approach is 
often necessary just to compute the inverse 
transform after a complicated substitution like 
BDF2. (Applying BDF2 to the kernel implied by 
Equation 16, a mere delay, results in a formula 
involving Hermite polynomials of high order!) 
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(b) 


Fig. 2. The response Z, of Eq. (21) in (a) free space, and (b) a 
Debye material, at different distances from a one-dimensional 
source. 


To demonstrate this computation, consider the 
Laplace domain kernel 





Ap 
1+ st 


Z(s) = exp| —sx [Po + , (25) 


which is related to the response of a Debye 
medium in one dimension. (In system engineering 
terms, a scalar function of the Laplace transform 
would be called a transfer function.) Fig. 2(a) 
shows the response of this medium for p,. = 1 and 
Ap = 0, and Fig. 2(b) shows the response of the 
medium for for p, = 1, Ap = 1, andt=1. (The 
computation was done with an FFT of M = 256 
points.) The dispersive effect of the Debye 
medium is clear from the lower peaks of the 
various curves, the damping of the oscillating 
approximation of the Green’s function, and the 
fact that the peak of the curves in the dispersive 
graph are to the right of those in the non-dispersive 
graph. 


D. A Runge-Kutta Based Approach 


The discussions of the previous subsections 
illustrate the CQ method of discretization and 
explain how it can be used to model complex 
kernels. The methods introduced so far have one 
important failing however: no A-stable multistep 
method can ever have accuracy greater than 
second order [25]. Fortunately, the remedy for this 
is also well known: High order solutions to stiff 
differential equations are provided by implicit 
Runge-Kutta (IRK) methods [23, 24, 26]. 


Most readers will recall the scheme known 
universally as “the” Runge-Kutta (RK) solution 
scheme [24, 26]. In this fourth-order scheme, the 
solution is advanced from a single previous time 
step to the next time step by using two 
successively computed approximations to the 
solution at the midpoint of the interval, and the 
combination of these approximations with the 
previous time step value to find the next 
approximation. This is a Runge-Kutta scheme 
because unlike multistep methods, it uses 
approximations internal to the interval to create the 
next approximation rather than the value at 
previous steps. It 1s an explicit method because 
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each approximation is created from the previous 
approximations successively without solving a set 
of equations. Thus, the approximations internal to 
the interval do not depend on the approximations 
on the right endpoint of the interval, and the two 
approximations at the interval midpoint are 
computed consecutively. In an implicit method, 
all of the approximations for a given time step are 
computed simultaneously by solving a set of 
equations. 


Even though “the” Runge-Kutta method is the 
most well known of the Runge-Kutta methods, 
there is a plethora of methods both explicit and 
implicit for solving ordinary differential equations 
[24]. In general, a p-stage RK method is a method 
to solve a general first order initial value problem 


d 
= = f(t) 


specified by a pxp matrix A = |a; jl, two p- 
vectors b = [b,;| and c = [c;], and a time step At. 
(The parameters A, b, and c are often collectively 
called the Butcher tableau of the RK method [24, 
26].) Given the approximation from a previous 
step Yn ~ y(nAt), the intermediate values Y,,; 
(i = 1,...,p) needed to find Yp+ı are found by 
solving the system of equations 
p 
Yni = Yn + at) aij f (tn + c,At, Ynj): 


j=l 


y(0) = yo (26) 


(27) 


Once the Y,; are known, the approximation to the 
solution at the next time step is computed using 
the formula 

p 


Yn+1 ~ Yn + at) bj f (tn $ cjAt, ar) (28) 
j=l 
Explicit methods have a strictly lower triangular A 
matrix with a vanishing diagonal. Implicit methods 
usually have a fully populated matrix. The only 
methods that are A-stable are IRK methods [24, 
26]. 

There are a host of ways of discovering IRK 
methods with different accuracy and stability 
properties, as well with different node locations. 
For our use here, the most important IRK methods 
are those known as the Radau-IJA methods. 
Radau integration rules are similar to the more 


familiar Gauss-Legendre integration rules, but they 
are constructed to force a single endpoint of an 
interval to be one of the nodes of integration. 
(This in turn implies they have slightly lower order 
than the corresponding Gauss-Legendre rule for 
the same number of integration points.) In Radau- 
IHA IRK methods, Radau integration rules are used 
to create an A-stable method in which ¥,, = Yn+1; 
that is, the last internal node is also the next 
function approximation. In terms of the Butcher 
tableau of the Radau-IIJA methods, this means that 
Cp = 1 and a,; = bj. Butcher tableaus for two 
such methods are shown in Fig. 3. 





























5 1 i 
12 12 = 
A» = 3 1 bz = g 
Ža Ae 1 
88 -7V6 296—-169V6 -2 +3vVv6 Z 
360 1800 225 4 a Q 
296+169V6  88-7V6  -—2- 3V6 
1800 360 225 
16=46 1644/6 1 = 
36 36 9 


Fig. 3. The Butcher tableaus of two- and three-stage Radau IIA 

methods. 

To apply IRK methods to the discretization of 
our integral equations, we need to find the 
appropriate substitution for s in terms of z. To this 
end, we apply Eqs. (27) and (28) to the solution of 
the canonical Eq. (14) in the Z transform domain. 
If we collect the values Y,,; into a p-vector Y,,, and 
then in the Z domain, we find from Eq. (28) 

SING nee 
y(z) = -b YC) (29) 
where the superscripted “T” indicates the 
transposition of a column vector into a row vector. 
Treating Eq. (27) similarly results in 


Y(z) = uy(z) + sAtAY(z). (30) 


Solving this equation for y(z), and substituting the 
result back into Eq. (29) yields [15] 

T 
i + A Y(z). 





F(z) = sAt | nb (31) 


Z — 
In view of the quantity multiplying Y(z) on the 
right hand side of this equation, we find the matrix 
substitution 
-1 
+ A 





(32) 


= At 


1 ub! 
z— 1 
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for s in the Laplace domain integral equation. 
This equation is a finite difference method in the 
frequency domain, but it is disguised: More than 
one time step is being computed at once, and the 
internal time step size is an irrational fraction of 
At. 


E. Dispersion Mitigation 


The CQ method as described up until this point 
has many attractive properties: It provides a 
temporal discretization that is easy to integrate 
spatially, it can be easily applied to situations with 
complicated Green’s functions, and it has an 
accuracy that is adjustable to yield very fast 
convergence with shrinking time step. 


Most of these benefits are directly related to the 
biggest drawback of the method. As shown in Fig. 
1, the method has a built-in numerical dispersion 
that affects detrimentally both the speed and 
accuracy of the result. Of course, this is primarily 
a concern in non-dispersive media, where the 
compromised accuracy is noticeable and the 
computational deceleration is easily noticed. 


One way to mitigate numerical dispersion and 
dissipation is to use a higher order time stepping 
method (or a smaller timestep!). This delays the 
onset of serious inaccuracy. Fortunately, in non- 
dispersive media, another fix presents itself 
immediately: Because the primary salubrious 
ramification of the fictional dispersion 1s the ease 
of integration in the near field, the dispersion and 
dissipation can be halted once this effect ceases to 
be so helpful. To see how this can be done, 
consider again the discretization of the Laplace 
domain kernel of Eq. (18) for time step At and 
different values of the delay t. Choose an integer 
time step threshold K after which dispersion will 
be halted. To discretize the delay, follow the 
standard CQ prescription if tT < KAt. Ift > KdAt, 
compute the integer 

M = |- - K| (33) 

At 
where [x| is the smallest integer greater than or 
equal to x. The discretization of the delay can 
now be computed by using the CQ recipe not with 
the actual delay T, but with the much shorter delay 


T! =T — Måt, (34) 


a number between (K — 1)At and KAt. Once the 
CQ discretization of t’ is computed, the result can 
be delayed M time steps to yield the final 
approximation. 


Fig. 4 demonstrates the effect of applying the 
BE method both with and without dispersion 
mitigation. Specifically, Fig. 4(a) shows the 
discretization of the delay kernel of Eq. (18) with a 
discretization size At = 1 for values of T ranging 
from 0.25 to 40. The graphs are stacked vertically, 
with the largest delay at the top. The effects of 
dispersion and dissipation are clearly visible in this 
picture, as the topmost line 1s almost completely 
flat. 
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Fig. 4. The temporal evolution of a delay discretized by (a) standard 
BE, and (b) BE with dispersion halted after K = 5 time steps. 
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Fig. 4(b) shows plots constructed with the same 
parameters, but with dispersion and dissipation 
halted at K = 5 time steps. The almost complete 
lack of change in the shape of the curve with 
increasing delay after a short preliminary phase is 
evident. 


Il]. NUMERICAL RESULTS 


In this section, we apply the techniques of the 
previous sections to scattering problems to 
demonstrate their accuracy and convergence. In 
each example, the x-polarized, z-directed incident 
wave is of the form 


pe exp (t= tp) teost 20 fot |, (35) 


where t = t — z/c. The frequency content of this 
wave is controlled by fọ, which sets the center 
frequency of the incident wave, and o, which sets 
its bandwidth around that frequency. In particular, 
we define o through a nominal bandwidth B 
(measured in Hz) with the formula 
3 
-nB 
This definition makes the power contained outside 
the frequency range [fp — B, fo + B| unobservable 
at double precision. The parameter t, is a delay 
incorporated so the maximum field reaches the 
origin at a positive time. Because fọ + B is 
considered the maximum frequency contained in 
the incident wave, results are related in terms of an 
oversampling factor y% defined by 
1 
Y = Oh + BAe’ 
0 

that is, w is the ratio of the sampling frequency 
used to the Nyquist frequency, assuming that the 
incident wave is actually bandlimited to its 
nominal maximum frequency. 


(36) 


(37) 


The first example examined here demonstrates 
the basic ideas for the scattering from a conducting 
sphere of diameter Im, excited by a wave with 
center frequency fp = 200 MHz, with B = 150 
MHz, and t, = 55 ns. The sphere was meshed by 
projecting an octahedron meshed into 128 
equilateral triangles onto the surface of the sphere. 
(The resulting triangles are of course perfectly 


spherical.) Graglia, Wilton, Peterson basis 
functions of first order (4.e., quadratic functions) 
were used for the modeling, resulting in 640 
unknowns. Simulations were done for various 
oversampling factors, and the bistatic radar cross 
section was computed for 41 frequencies equally 
spaced between 100 MHz and 300 MHz inclusive, 
and 91 values of azimuth between 0 and 180 
degrees. 
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Fig. 5. Bistatic RCS ofa 1m diameter conducting sphere computed 
in the time domain using four different methods at (a) 100 MHz and 
(b) 300 MEZ. 

Fig. 5 shows the scattering from the sphere using 
four different temporal discretizations. The BE 
and BDF2 discretizations were computed with a 
time step size of At = 200 ps, the two-stage Radau 
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method was computed with At = 400 ps, and the 
three stage Radau result was computed with t = 
500 ps. (Of course, it must be remembered that 
the n-stage Radau method has n samples per time 
step.) The time step values represent an 
oversampling factor just above 7 for BE, BDF2, 
and the two-stage RadaullA formula, and slightly 
more for the three-stage RadaulIA formula. 
Convergence of all four methods to a frequency 
domain method of moments (MoM) result is 
shown in Fig. 6. 
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Fig. 6. Convergence of bistatic RCS (measured in dB and 
computed in the one norm) to the frequency domain result for four 
different CQ-based methods for a PEC sphere. 


To demonstrate the technique of dispersion 
mitigation, a single triangular patch touching the 
pole on the negative z-axis was removed from the 
sphere that was excited by the same wave. The 
coupling between the incident wave and this small 
hole was not large, but allowed for enough 
oscillation to see an effect. Fig. 7 shows the 
magnitude of the current density at the centroid of 
a patch computed both with the straightforward 
CQ method, and also dispersion limited with 
K = 4. Fig. 8 shows the convergence of the 
solution in both cases. These figures show that 
while the dispersion limitation scheme reduces the 
fictitious current that oscillates once the simulation 
falls below its accuracy threshold, this does not 
much alter the overall accuracy or convergence 
rate. Of course, the dispersion-limited simulation 
delivers results much faster than the unlimited 
simulation. 


Next, the same spherical geometry was used to 
demonstrate the ability of the CQ method to 
simulate physical dispersion. For this example, 
the dielectric constant of the sphere was assumed 
to have a Debye pole and frequency domain 
behavior given by 

1.75 


1+5x10-1%s 
The sphere was again excited by the same incident 
wave. Fig. 9 shows scattering results for the same 
strategies and time steps as Fig. 5. Fig. 10 shows 
the convergence of the scheme. 


e(s) = 2.25 + (38) 
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Fig. 7. Currents returned by dispersion limited and dispersion 

unlimited BDF2 simulations. 
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Fig. 8. Convergence of bistatic RCS of a sphere to a frequency 
domain result for time domain simulations performed with and 
without dispersion limitation for a sphere with a triangular hole. 
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Fig. 9. Bistatic RCS of a 1m diameter penetrable sphere composed 
of Debye material computed in the time domain using four different 
methods at (a) 100 MHz and (b) 300 MHz. 


As a final example, we examine the scattering 
from a highly resonant structure to demonstrate the 
acceleration possible using the dispersion 
mitigation method of Section II-E. The structure 
in question is a rectangular prism of dimensions 
0.5mx0.5mx1.0m. The square sides are parallel to 
the x-y plane, and the entire structure sits in the 
octant in which all of the coordinates are positive 
with one corner at the origin. On the bottom 
square, a 0.1m wide slot is cut in the prism, 
centered in the y-direction and running the length 
of the prism in the x-direction. 


This structure was excited with the Gaussian 
plane wave of Eq. (35) with fọ = 300 MHz, 
B = 200 MHz, and t,=55 ns, and the 
backscattering from the structure was computed at 


91 frequencies evenly spaced between 150 and 
450 MHz, inclusive. 
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Fig. 10. Convergence of different simulation schemes of bistatic 
RCS of a Debye dielectric sphere. 


Figure 11 shows the monostatic RCS of this 
scatterer, illuminated by a wave aimed at its slot, 
and polarized parallel to it. The results follow the 
frequency domain MoM well (so long as the 
method employed is better than first-order 
accurate) even given the obviously resonant 
behavior on display as shown in Figure 12. 
Figures 13 and 14 compare the setup and marching 
time for the dispersion arrested and standard CQ 
schemes. The dispersion-arrested scheme is faster 
to setup and march than the pure CQ scheme. 
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Fig. 11. Monostatic RCS of a slotted, conducting prism computed 
by several methods. 
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IV. CONCLUSIONS 


CQ methods are an interesting alternative to TG 
methods for the temporal discretization of time 
domain integral equations. They easily achieve 
very high accuracy, ensure stability over a wide 
range of parameters, and can be implemented for 
even the most complex Green’s functions. 
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(b) 
Fig. 12. Convergence of four time domain methods to MoM results 
for the slotted rectangular prism with (a) K = 4, and (b) K = œ. 


The primary drawback of CQ methods is the 
manner in which they introduce spurious 
numerical dispersion, which affects both their 
accuracy and their efficiency. This study 
introduced a method to mitigate the effects of this 
dispersion by halting the evolution of the radiated 
wave after a fixed number of time steps. While 
this greatly shortened run times, it failed to have 
much effect on accuracy. 


Future work will concentrate on combinations of 
this method with TG methods so that each method 
may be used where its properties are most 
advantageous. 
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(b) 
Fig. 13. Setup and marching time for two different difference 
formula based time stepping methods simulating the scattering from 
a slotted rectangular prism. 


REFERENCES 


[1] C. L. Bennett and W. L. Weeks, “Transient scattering 
from conducting cylinders,” IEEE Transactions on 
Antennas and Propagation, 18, pp. 627-633 (1970). 


[2] A. J. Poggio and E. K. Miller, “Integral equation 
solutions of three dimensional scattering problems,” in 
Computer Techniques for Electromagnetics, R. Mittra, 
Ed. Oxford, UK: Pergammon Press, 1973. 


12 


Forum for Electromagnetic Research Methods and Application Technologies (FERMAT) 


10° 





— RadaullIA-2 Setup 

-~=--RadaullA-2 March 

—— RadaullA-2 Setup (K= 4) 
RadaullA-2 March (K= 4) 


Time (s) 
= 

















l 10 100 


(a) 
10° 
—=— RadaullA-3 Setup 


--©--RadaulJA-3 March 
—<— RadaullA-3 Setup (K= 4) 









--X-- RadaullA-3 March (K= 4) 





PED 
2 w oe 
= ® 
= Pa 














1 10 100 
Y 


(b) 
Fig. 14. Setup and marching time for two different Runge-Kutta 
based time stepping methods simulating the scattering from a 
slotted rectangular prism. 


[3] S. M. Rao and D. R. Wilton, “Transient scattering by 
conducting surfaces of arbitrary shape,” JEEE 


Transactions on Antennas and Propagation, 39, 56-61 
(1991). 


[4] B. P. Rynne and P. D. Smith, “Stability of time 
marching algorithms for the electric field integral 
equation,” Journal of Electromagnetic Waves and 
Applications, 4, 1181-1205 (1990). 

[5] M. J. Bluck and S. P. Walker, “Time-domain BIE 
analysis of large three dimensional electromagnetic 
scattering problems,” IEEE Transactions on Antennas 
and Propagation, 45, pp. 894-901 (1997). 

[6] G. Manara, A. Monorchio, and R. Reggiannini, “A 
space-time discretization criterion for a stable time- 
marching solution of the electric field integral equation,” 


IEEE Transactions on Antennas and Propagation, 45, 
pp. 527-532 (1997). 


[7] S. M. Rao and T. K. Sarkar, “Implicit solution of time- 
domain integral equations for arbitrarily shaped 
dielectric bodies,” Microwave and Optical Technology 
Letters, 21, 201-205 (1999). 


[8] D. S. Weile, G. Pisharody, N.-W. Chen, B. Shanker, and 
E. Michielssen, “A novel scheme for the solution of 
time-domain integral equations,” IEEE Transactions on 
Antennas and Propagation, 52, 283-295 (2004). 


[9] T. Abboud, J.-C. Nédélec, and J. Volakis, “Stable 
solutions of the retarded potential time domain integral 
equations,” Applied Computational Electromagnetics 
Conference (ACES) 2001. 


[10]T. Ha-Duong, B. Ludwig, and I. Terasse, “A Galerkin 
BEM for transient acoustic scattering by an absorbing 
obstacle,” International Journal for Numerical Methods 
in Engineering, 57, 1845-82 (2003). 

[11]B. H. Jung, Z. Ji, T. K. Sarkar, M. Salazar-Palma, andM. 
Yuan, “A comparison of marching-on in time method 
with marching-on in degree method for the TDIE 
solver,” Progress in Electromagnetic Research (PIER), 
70, 281-296 (2007). 


[12] A. J. Pray, Y. Beghein, N. V. Nair, K. Cools, H. Bagcı, 
B. Shanker, “A Higher order space-time Galerkin 
scheme for time domain integral equations,” IEEE 


Transactions on Antennas and Propagation, 6183-6191 
(2014). 


[13]C. Lubich, “Convolution quadrature and discretized 
operational calculus I,” Numerische Mathematik, 52, 
129-145 (1998). 


[14] C. Lubich, “Convolution quadrature and discretized 
operational calculus I,” Numerische Mathematik, 52, 
413-425 (1998). 

[15] C. Lubich and A. Ostermann, “Runge-Kutta methods 
for parabolic equations and convolution quadrature,” 
Mathematics of Computation, 60, 105-131 (1993). 

[16]W. Hackbusch, W. Kress, and S. Sauter, Sparse 
convolution quadrature for time domain boundary 
integral formulations of the wave equation, Tech. Rep. 
116, Max Planck Institute for Mathematics in the 
Sciences, Leipzig, Germany (2005). 

[17]B. P. Lathi, Linear Signals and Systems, 2" Revised 
Edition, Oxford University Press (2009). 

[18] A. V. Oppenheim and R. W. Schaeffer, Discrete Time 
Signal Processing, 3™ Edition, Prentice-Hall (2009). 
[19] W. L. Brogan, Modern Control Theory, 3’ Edition, 

Prentice Hall (1991). 

[20]R. F. Harrington, Field Computation by Moment 
Methods, Robert E. Krieger Publishing Company 
(1968). 

[21] A. F. Peterson, S. Ray, and R. Mittra, Computational 
Methods for Electromagnetics, IEEE Press (1998). 

[22] X. Wang, R. A. Wildman, D. S. Weile, and P. Monk, “A 
finite difference delay modeling approach to the 


13 


Forum for Electromagnetic Research Methods and Application Technologies (FERMAT) 


discretization of the time domain integral equations of 
electromagnetics,” IEEE Transcations on Antennas and 
Propagation, 56, 2442-2452 (2008). 

[23]X. Wang and D. S. Weile, “Implicit Runge-Kutta 
methods for the discretization of time domain integral 
equations,” IEEE Transcations on Antennas and 
Propagation, 59, 4651-4663 (2011). 

[24]K. E. Atkinson, W. Han, and D. E. Stewart, Numerical 
Solution of Ordinary Differential Equations, John Wiley 
and Sons (2009). 


[25]G. Dahlquist, “Convergence and stability in the 
numerical integration of ordinary differential equations, 
Mathematica Scandinavica, 4, 33-53 (1956). 


[26]J. C. Butcher, Numerical Methods for Ordinary 


Differential Equations, 2" Edition, John Wiley and Sons 
(2008). 


Jielin Li received the B.S. 
degree in electronic information 
science and technology from 
Wuhan University, Hubei, China 
in 2008. Currently he is a 
graduate research assistant with 
the Department of Electrical and 
Computer Engineering, University of Delaware. 
His current research interests are in computational 
electromagnetics, especially time domain integral 
equations. At present, he is working on the fast 
method to solve the time domain integral 
equations. 








MPH Peter Monk received the B.A. 
WZ degree from Cambridge University, 
U.K., in 1978 and the Ph.D. degree 
from Rutgers University, USA in 
1983. Since 1982 he has been on 
the faculty of the Department of 
Mathematical Sciences at the University of 
Delaware in the USA. He was appointed Unidel 
Professor in 2000. His research interests include 
applications of finite element methods to the 
simulation of solar-voltaic structures, Trefftz 
methods and inverse electromagnetic scattering. 
He is the author of roughly 130 papers and a book 
on “Finite Element Methods for Maxwell’s 
Equations.” As a result of delivering a series of 
lectures sponsored by CBMS-NSF he has also 
published a research monograph entitled “The 
Linear Sampling Methods in _ _ Inverse 
Electromagnetic Scattering” coauthored with F. 








Cakoni and D. Colton. 


Daniel S. Weile obtained his 
B.S.E.E and his B.S. (in 
Mathematics) at the University of 
Maryland at College Park and 1994, 
and M.S. and Ph.D. in Electrical 
Engineering at the University of 
Illinois at Urbana-Champaign in 1995 and 1999, 
respectively. | Currently, he is an Associate 
Professor of Electrical Engineering at the 
University of Delaware. In 1994, he worked at the 
Institute for Plasma Research developing 
interactive software for the design of depressed 
collectors for gyrotron beams. As a research 
assistant and Visiting Assistant Professor at the 
University of Illinois, Dr. Weile worked on the 
efficient design of electromagnetic devices using 
stochastic optimization techniques, and fast time- 
domain integral equation methods for the solution 
of scattering problems. His current research 
interests include computational electromagnetics 
(especially time-domain integral equations), 
periodic structures, and the use of evolutionary 
optimization in electromagnetic design. Dr. Weile 
is the recipient of an NSF CAREER Award and 
and an ONR Young Investigator Award. He is a 
member of Eta Kappa Nu, Tau Beta Pi, Phi Beta 
Kappa, and URSI Commission B. 


14 


